library(CellChat)
library(patchwork)
library(tidyverse)
library(Seurat)
library(tidyseurat)
library(ggpubr)
source('00_util_scripts/mod_seurat.R')

crc4 <- read_rds('CRC-I/data/crc_merge4_immune.rds')

augur_crc4 <- crc4 |>
  mutate(tt = genotype == 'TT') |>
  calculate_auc(label_col = 'tt', cell_type_col = 'manual_main')


obcc.tt.crc4 <- crc4 |>
  filter(genotype == 'TT' & !str_detect(manual_main, 'Dividing')) |>
  createCellChat(group.by = 'manual_main')

obcc.ctrl.crc4 <- crc4 |>
  filter(genotype != 'TT' & !str_detect(manual_main, 'Dividing')) |>
  createCellChat(group.by = 'manual_main')

# set ppi db
showDatabaseCategory(CellChatDB.human)

obcc.tt.crc4@DB <- obcc.ctrl.crc4@DB <- subsetDB(CellChatDB.human)

# infer cell-cell network ------------
obcc.tt.crc4 <- obcc.tt.crc4 |>
  infer_cellschat()

obcc.ctrl.crc4 <- obcc.ctrl.crc4 |>
  infer_cellschat()

obcc_list <- list(
  'II+IT' = obcc.ctrl.crc4,
  'TT' = obcc.tt.crc4
)

obcc <- obcc_list |>
  mergeCellChat(add.names = c('II+IT', 'TT'))

# compare total interaction number & strength --------
g1 <- obcc |>
  compareInteractions(group = c('II+IT', 'TT'),
                      color.use = c('blue','red'))

g2 <- obcc |>
  compareInteractions(group = c('II+IT', 'TT'),
                      color.use = c('blue','red'), measure = 'weight')

g1 + g2 & theme_pubr(legend = 'none')

# diff among net ---------
## chord plot
par(mfrow = 1:2, xpd = T)
obcc |> netVisual_diffInteraction(weight.scale = T)
obcc |> netVisual_diffInteraction(weight.scale = T, measure = 'weight')

## heatmap
obcc |> netVisual_heatmap()
obcc |> netVisual_heatmap(measure = 'weight',
                          font.size = 12,
                          font.size.title = 12)

## incoming & outgoing --------------
num.link <- obcc_list |>
  sapply(\(x){
    y <- x@net$count
    rowSums(y) + colSums(y) - diag(y)
  })

weight.min.max <- c(min(num.link), max(num.link))  

gg <- list()

for (i in 1:2) {
  gg[[i]] <- obcc_list[[i]] |>
    netAnalysis_signalingRole_scatter(title = names(obcc_list)[i],
                                      weight.MinMax = weight.min.max)
}

wrap_plots(gg)

### identify specific signal in one cell type
obcc |>
  netAnalysis_signalingChanges_scatter(idents = 'B',font.size = 12,
                                       color.use = c('grey','blue','red')) +
  theme_pubr(legend = 'right')

### information flow of pathways ---------
gg1 <- obcc |>
  rankNet(do.stat = T, color.use = c('blue','red'),font.size = 10) +
  NoLegend()

gg2 <- obcc |>
  rankNet(do.stat = T, color.use = c('blue','red'), font.size = 10, stacked = T)

gg1 + gg2

### outgoing pattern ----------
pathway.union <- union(obcc_list[[1]]@netP$pathways,
                       obcc_list[[2]]@netP$pathways)

ht1 <- obcc_list[[1]] |>
  netAnalysis_signalingRole_heatmap(pattern = "outgoing",
                                    signaling = pathway.union,
                                    title = names(obcc_list)[1],
                                    width = 5, height = 6)

ht2 <- obcc_list[[2]] |>
  netAnalysis_signalingRole_heatmap(pattern = "outgoing",
                                    signaling = pathway.union,
                                    title = names(obcc_list)[2],
                                    width = 5, height = 6)

ComplexHeatmap::draw(ht2 + ht1, ht_gap = unit(1, "cm"))

### differential ligand pairs -----------
obcc |>
  netVisual_bubble(angle.x = 45,
                   sources.use = 5,
                   targets.use = 2:3,
                   comparison = 1:2,
                   color.text = c('blue','red'),
                   max.dataset = 2,
                   remove.isolate = T,
                   title.name = 'Increased signaling in TT')

obcc |>
  netVisual_bubble(angle.x = 45,
                   sources.use = 5,
                   targets.use = 2:3,
                   comparison = 1:2,
                   color.text = c('blue','red'),
                   max.dataset = 1,
                   remove.isolate = T,font.size = 12,
                   title.name = 'Decreased signaling in TT')

#### use DEA method ---------
obcc <- obcc |>
  identifyOverExpressedGenes(group.dataset = "datasets",
                             pos.dataset = 'TT',
                             features.name = 'TT.merged',
                             only.pos = FALSE, thresh.pc = 0.1,
                             thresh.fc = 0.05,thresh.p = 0.05,
                             group.DE.combined = T)

net <- obcc |>
  netMappingDEG(features.name = 'TT.merged',
                variable.all = T)

net.up <- obcc |>
  subsetCommunication(net = net, datasets = "TT",
                      ligand.logFC = 0.05, receptor.logFC = NULL)

net.down <- obcc |>
  subsetCommunication(net = net, datasets = "II+IT",
                      ligand.logFC = -0.05, receptor.logFC = NULL)

pairLR.use.up = net.up[, "interaction_name", drop = F]

gg1 <- obcc |>
  netVisual_bubble(pairLR.use = pairLR.use.up,
                   sources.use = 4,
                   targets.use = c(5:11),
                   comparison = c(1, 2),
                   angle.x = 90,
                   remove.isolate = T,
                   title.name = paste0("Up-regulated signaling in ", names(obcc_list)[2]))

pairLR.use.down = net.down[, "interaction_name", drop = F]

gg2 <- obcc |>
  netVisual_bubble(pairLR.use = pairLR.use.down,
                   sources.use = 4,
                   targets.use = c(5:11),
                   comparison = c(1, 2),
                   angle.x = 90,
                   remove.isolate = T,
                   title.name = paste0("Down-regulated signaling in ", names(obcc_list)[2]))

gg1 + gg2

# Chord diagram
par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_gene(object.list[[2]], sources.use = 4, targets.use = c(5:11), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 4, targets.use = c(5:11), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))
#> You may try the function `netVisual_chord_cell` for visualizing individual signaling pathway


